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The Boussinesq equations for Rayleigh-Benard convection are simulated for a cylindrical con- 
tainer with an aspect ratio near 1.5. The transition from an axisymmetric stationary flow to time- 
dependent flows is studied using nonlinear simulations, linear stability analysis and bifurcation 
theory. At a Rayleigh number near 25 000, the axisymmetric flow becomes unstable to standing 
or travelling azimuthal waves. The standing waves are slightly unstable to travelling waves. This 
scenario is identified as a Hopf bifurcation in a system with 0(2) symmetry. 



1. Introduction 

Rayleigh-Benard instability in a fluid layer heated from below in the presence of gravity is 
the classic prototype of pattern formation. A new chapter in its investigation began with the in- 
crease of computer performance that made feasible three-dimensional, nonlinear, high-resolution 
simulations of the Boussinesq equations governing this system. 

We are interested in a fluid layer confined in a vertical cylinder whose upper and lower bound- 
ing surfaces are maintained at a temperature difference measured by the Rayleigh number. The 
conductive solution for this system is a motionless state with a uniform vertical temperature 
gradient. This solution is stable up to a critical Rayleigh number Ra c , whose value depends on 
the aspect ratio F = radius/height. Above Ra c , convective motions appear and form various roll 
structures. 

A summary covering the develop ments since the mid 1980 s for conv ective systems with large 
aspect ratio (T 3> 1) can be found in lBodenschatz. Pesch & Ahlersl lEoOO). In such domains a rich 



variety of patterns was reported: "P an Am" patterns (arche s with several ce ntres of curvature, see 
Ahlers. Cannell & Steinberg! 1985h. straight parallel rolls dCroquettel[l 989: Croquette, Le Gal & 
Pocheau ll986h~ concentric rolls (targets, see Koschmieder & Pallaslll974l: Croquette, Mory & 



Schos seler 19831) , one- and several- armed rotating sp irals dPlapp. Egolf. Bodenschatz & Pesch 
1 1998h . target s with dislocated centre dCroauettell 19891). hexagonal cells (Ciliberto, Pampalon i & 
Perez-Garcfa 1988b and spiral-defect chaos dMorris. Bodenschatz. Cannell & Ahlers 1993 ). A 



large ov erview on convective phenomena observed experimentally before this time can also be 



found in lKoschmied er ( 1993). 



We focus here on cylinders with moderate aspect ratio r ~ 1 . The flow structure then depends 
strongly on system geometry. Fo r this regime, the stabilit y of the conductive state was well es- 
tablish ed in the 1970s-1980s bvlCharlson & Sanil dl970h. IStork & Mullerl dl975l) and Buell & 
C atton d 19831) . Critical Rayleigh numbers Ra r are about 2000 for ~ > 1, increasing steeply for 



lower r and decreasing asymptotically towards Ra c = 1708 for r — > °°. Charlson & Sani ( 1970h 
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estimated by a numerical variational technique the onset of axisymmetric convection in cylin- 
ders of aspect ratios between 0.5 and 8, with insulating and conducting sidewalls. They found 
the critical Rayleigh numbers (Ra c = 2545 for r = 1, decre asing for higher F) and the corre- 



sponding number of rolls. They then generalised this analysis dCharlson & Sani|[l971 ), including 



non-axisymmetric mo des and predicting Ra c and corresponding critical azimuthal wavenumbers. 



Stork & Miillerl ( 1975 ) observed experimentally convective patterns in annuli and cylinders of 
aspect ratio 0.7 < F < 3.2, varying the sidewall insulation. Their critical Ray l eigh n umbers were 
in good agreement with those predicted by Charlson and Sani. iRosenblatl dl982l) investigated 
convective instabilities numerically for free-slip boundary conditions, using a severely truncated 
expansion in a small number of eigenmodes. He described n on-axisymmetric motio ns existing 
just above onset for aspect ratios between 0.5 and 2.0. Finally. IBuell & Cattonl (119831) described 
how the onset of convection is influenced by the ratio of the fluid conductivity to that of the wall, 
by performing linear analysis for the aspect ratio range < F < 4. They determined the critical 
Rayleigh number and azimuthal wavenumber as a function of both aspect ratio and sidewall con- 
ductivity, thus completing the results of the previous investigations, which considered either per- 
fectly insulating or perfectly conducting walls. These results were confirmed bv Marques, Net, 
Massaguer & Mercader dl993l) . The flow succeeding the condu ctive state is three-d imensional 



over large ranges of aspect ratios, contrary to the expectations of lKoschmiedei ( 19931) . 



The stability of the first convective state, depending on both aspect ratio and Prandtl number, 
has been investigated mainlv for situations in which the primary flow is axisvmmetric. Charlson 
& Sani dl975h attempted to predict numerically the stability of the primary axisymmetric flow, 



but the resolution available at that time was inadequate to the task. Miiller, Neumann & Weber 



(1984) investigated convective flows experimentally and theoreticall y. They observed axi sym 



metric flows for F =1 and non-axisymmetric flows for 0.1 < F < 0.5. Hardin & Sani ( 1993|) cal- 
culated weakly nonlinear solutions to the Boussinesq equations for several moderate and small 
aspect ratios. They found a bifurcation from the axisymmetric state towards a mode with azimu- 
thal wavenumber m = 2 for F = 1, Pr = 6.7 and Ra C 2 — 2430. 

The most complete numerical study of secondary convective i nstabi lities for moderate aspect 
ratio cylinders was performed by lWanschura. Kuhlmann & Rathl (1996). For cylinders with insu- 
lating sidewalls and 0.9 < F < 1 .57, the primary bifurcation to convection occurs at Ra c i=s 2000 
and leads to an axisymmetric flow whose stability was investigated for Prandtl numbers 0.02 
and 1 . Wanschura et al. predicted the succeeding flows to be steady, except over a narrow aspect 
ratio range 1 .45 < F < 1 .57 at Pr = 1, where they found oscillatory instabilities at Ra c 2 « 25000 
towards flows with azimuthal wavenumbers m = 3 and m = 4. The primary aim of this paper is 
to provide a more detailed description of these bifurcations. 



Touihri, Ben Hadid & Henry d 19991) numerically investigated the stability of the conductive 



state for aspect ratios F = 0.5 and r = 1. They described the main critical modes and established 
a diagram of primary bifurcations, including unstable branches. They also found a secondary 
bifurcation point Ra c 2, at which the axisymmetric flow becomes unstable towards a two-roll flow 
and calculated Ra c a for r = 1 and < Pr < I. 

An interesting experimental study was carried out by iHof. Lucas & Mu llin dl999h . Varying 



the Rayleigh number through different sequences of values, for fixed parameters F = 2.0 and 
Pr = 6.7, they obtained several different stable patterns for the same final Rayleigh number. 
They also reported a transition from an axisymmetric steady state towards azimuthal waves. Our 
numerical simulations of this phenomenon are the subject of a separat e investigation. 

Mor e recently con vective patterns were numerically investigated by Rudige r & Feudel (2000) 



and bv lLeong] ([2002). Rudiger and Feudel found stability ranges for multi-roll patterns, targets 
and spirals for r = 4, Pr = I. Leong observed several steady convective patterns for aspect ratios 
2 and 4 and Prandtl number Pr — 7, all of which were stable in the range 6250 < Ra < 37 500, 
and calculated the heat transfer for each pattern. 
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FIGURE 1 . Geometry and coordinate system. 



Convective systems often display oscillatory behavior. In binary fluid or rotating convection, 
the primary bifurcation is usually to periodic states, while in Rayleigh-Benard convection, perio- 
dic behavior occurs as a secondary bifurcati on. The oscillatory and ske w- varicose instabi l ities o f 
long straight parallel rolls calculated in, e.g. lClever & Bussd (119741) and lBusse & Clever ( 1979), 
are manifested as travelling waves along rolls an d as periodic defect nucleation dCroquettel 19*891 : 
Croquette et al. l ll986tlRiidiger & Feudelll200oT) : rotating spirals were observed by the same in- 



vestigators; and radially propagating patterns of concentric rolls were observed by Tuckerman 
& Barkley il988j). However, none of these manifestations of oscillatory behavior resemble the 
azimuthal waves we describe in this study. 

Competition between standing and rotating azimuthal waves has been extensively studied in 
thermocapillary convection, driven by surface-tension gradients. For example, competition be- 
tween rotati ng and standing wave s is observed on the upper free surface of an open cylindrical 
container by ISim_& Zebib (2002 |) and in the midplan e of a cylindrical liquid bridge with free 
outer surface bv lLevpoldt. Kuhlmann & Rath ( 2000l) . both of aspect ratio 1. These azimuthal 
waves are very similar to those we describe in this study; however, such flows are uncommon in 
the Rayleigh-Benard (buoyancy-driven) convection literature. 

We wished to study in detail the time-periodic non-axisymm etric states in cylindric al Ray- 
leigh-Benard convection resulting from the bifurcation found bv lWanschura et al. (119961) . Hence 
we have simulated numerically the loss of stability of the first convective axisymmetric solution 
undergoing an oscillatory bifurcation for 1 .45 < Y < 1 .57 and Pr = 1 . In this paper we describe 
the results of nonlinear simulations and linear stability analysis, which identify the scenario in 
terms of bifurcation theory in systems with symmetries. 

In addition to obtaining results particular to cylindrical Rayleigh-Benard convection with 
these parameter combinations, our purpose is to demonstrate how numerical and theoretical 
techniques can be combined in order to obtain a complete bifurcation-theoretic understanding 
of the oscillatory states produced by this secondary bifurcation. Such an approach can be applied 
to analyse transitions in a wide vari ety of other physical systems, ranging from flows driven 
by differentia lly rotating boundaries (|Nore. T uckerman^ Daube & Xin 20031) to Bose-Einstein 
condensation (IHuepe. Tuckerman. Metens & Brachetll2003h . 



2. Method 

2.1. Governing equations 

We consider a fluid confined in a cylinder of depth d and radius R (figure [TJ. The aspect ratio 
is defined as Y = R/d. The fluid has kinematic viscosity v, density p, thermal diffusivity K and 
thermal expansion coefficient (at constant pressure) y. The top and bottom temperatures of the 
cylinder are kept constant, at 7b — AT/ 2 and 7b + AT/ 2, respectively, leading to the linear con- 
ductive temperature profile T(z) = To — zAT/d. The lateral walls are insulating. The Rayleigh 
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number Ra and the Prandtl number Pr are defined by 

(2.1a) 

KV 

Pr=-. (2.1b) 

K 

Using the units d 2 /K, d, K/d and VK/ygd 3 for time, distance, velocity and temperature, we define 
u and h to be the nondimensionalised velocity and deviation of the temperature from the basic 
vertical profile, respectively. We obtain the Boussinesq equations governing the system: 

Pr 1 (d,u + (u ■ V) u) = -Vp + Au + he z (2.2a) 
d,h + (u-V)h = Rau z + Ah (2.2b) 
V-u = 0. (2.2c) 

The boundary conditions for velocity are no-slip and no-penetration 

u = for r = T or z = ±l/2. (2.3) 

Since the horizontal plates are assumed to be perfectly conducting (Dirichlet condition for h) and 
the vertical walls are insulating (Neumann condition), the boundary conditions for the tempera- 
ture are 

h = for z = ±l/2, (2.4a) 
dh 

— = for r = T. (2.46) 
dr 

2.2. Symmetries 

Symmetries play an important role in the possible transitions undergone by this system. The 
Boussinesq equations ( 12.2b with boundary conditions ( 12.3b — ( f2~4T > have reflection symmetry in 
the vertical direction z, and rotational and reflection symmetry in the azimuthal direction 0. The 
reflection symmetry in z is broken by the first bifurcation to a convective state. If the first con- 
vective state consists of axisymmetric convective rolls, then its remaining symmetries are re- 
flection and rotation in 8, which together comprise the symmetry group 0(2). Bifurcations in 

large number 
van Gils & 

Mallet-Paret dl986l) : iKuznetsovl dl998l) : iCoullet & Ioossl dl990t> . We give a brief summary of 
their results. 

First, the critical eigenvector may be axisymmetric. This case may be further subdivided ac- 
cording to whether the eigenvector is reflection-symmetric or antisymmetric in and whether 
the eigenvalue is real or complex. A reflection-symmetric eigenvec tor can lead to a target pat- 



the presence of 0(2 ) symmetry w ere studied and classifi e d duri n g the 1980s by a la 
of researchers. e.g.lBaiail dl982h:rGolubitskv & Stewart] dl985t): Knobloch (1986); 



tern of radially propagating rolls, e.g. lTuckerman & Barklev ( 1988 ). The breaking of reflection 



symmetry is associated with azimuthal flow. 

Secondly, the critical eigenvector may be non-axisymmetric. If the critical eigenvalue is real, 
then the resulting bifurcation is a circle pitchfork, leading to a "circle" of steady states para- 
metrised by phase. Each steady state is reflection symmetric in 8 (about some value 0o). If 
reflection symmetry is broken by a subsequent bifurcation, the scenario is that of a drift pitch- 
fork, leading to slow motion ("drift") along the circl e. A complex eig e nvalue corresponding 



to a non-axisymmetric eigenvector, like that found by IWanschura et al. (1996) for parameters 



1.45 < r < 1.57, Pr = 1, Ra > 23000, leads to a Hopf bifurcation which engenders three non- 
linear branches: standing waves, counterclockwise travelling waves, and clockwise travelling 
waves. The standing waves are reflection-symmetric in 8 (again about some value 8o), while the 
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travelling waves break this symmetry. Our aim is to determine which of these types of waves is 
realised by our physical system. 

2.3. Numerical integration 



We integrated the equations by a classical pseudospectral method ([Gottlieb & Orszagll 19771) . in 
which each scalar / of the fields u and h is represented using Chebyshev polynomials in the 
radial and vertical direction and Fourier series in the azimuthal direction 

N r .N-.N e 

f(r,z,B,t) = £ f jkm (t)Cj(r/r)C k (2z)e im6 + c.c, (2.5) 

j,k : m=0 

where the permit ted combina t ions o f (j,m) are restricted by the parity and regularity condi- 
tions described in iTuckerman (1989) for u r , we, u z and h. The nonlinear (advective) terms were 



calculated in physical space and integrated via the Adams-Bashforth formula, while the linear 
(diffusive) terms were calculated in spectral space and integrated via the Crank-Nicolson for- 
mula. An influence matrix method was used to impose incompressibility ( Tuc kermanll 1 9891) . A 
resolution of N r + 1 = 36, 2(Nq + 1) = 80, N z + 1 = 18 gridpoints or modes was found to be 
sufficient for nonlinear simulations. All computations were performed on the NEC SX-5 vector 
supercomputer, with time step 2 x 10~ 4 or 4 x 10~ 4 , depending on Ra, with CPU time per time 
step per grid point of 10~ 6 . 

2.4. Linear stability analysis 

An important additional element in understanding the phenomena undergone by the system is 
linear stability analysis. The procedure, which we summarise below, is described in more detail in 



mm 

Mamun & Tuckermanl l 19951) : ITuckerman & Barkle y (2000) and references therein. We linearise 



the equations about a steady state (U,H): 

Pr 1 (3 f u + (U ■ V) u + (u • V) U) = - Vp + Au + he z (2.6a) 
d,h + (U • V) A + (u ■ V)H = Ra u z + Ah (2.6b) 
V-u = 0. (2.6c) 

Equations (12.6b with boundary conditions (12.3b -( f2~4l are then integrated in time in the same way 
as the nonlinear equations ( 12.21 . We abbreviate the linear evolution problem ( 12.61 by 



(2.7) 



Temporal integration is equivalent to carrying out the power method on the approximate expo- 
nential operator, since 

' ;)(o. (2.8) 

In order to extract the leading real or complex eigenvalues (those of largest real part) and 
corresponding eigenvectors, we postprocess the results of integrating ( 12.6b as follows. A small 
number of fields 

^) (0) '(") (r) '(") (2r) '-'(") ((/ "~ 1)r) (2 ' 9) 

are calculated, by carrying out T /At linearised timesteps. The Krylov space corresponding to 
initial vector (u,h) T and matrix e LT is the ^-dimensional linear subspace consisting of all linear 
combinations of vectors in (12.9b . These vectors are orthonormalised to one another to generate a 
set of vectors Vi,V2,V3, ... vk which form a basis for the Krylov space. The action of the operator 
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on the Krylov space is represented by a small (K x K) matrix M whose elements are 

M jk = (vj,e LT v k ). (2.10) 

The small matrix M can be directly diagonalised. Its eigenvalues X approximate a small number 
of the eigenvalues of the large matrix e LT : this is the essence of Arnoldi's method. The procedure 
of generating the Krylov space via repeated action of e LT selects preferentially the K dominant 
values (those of largest magnitude) of e LT , i.e. the K leading eigenvalues (those of largest real 
part) of L. 

The eigenvectors of M prescribe coefficients of the vectors Vj which can be combined to 
form approximate eigenvectors § of e LT . The accuracy of these approximate eigenpairs (A,,<|)) 
is measured by the residue ||e Lr (j) — A-(j)|| in the case of real eigenvalues or by the residues 
| |e L7 >* - (X R <f - |, | |e L7 y - (A-V + A/(|) R )|| in the case of complex eigenvalues. If the 
desired eigenvalues have sufficiently small residues, they are accepted; otherwise we continue 
integration of ( 12.6b . replacing ( 12.91 ) by 

(;)(r) l (j)(2r) 1 (;)(3r) 1 ...,(j)(zr) (2.11) 

and so on, until the residue is below the acceptance criterion. 

After integrating the axisymmetric version of the nonlinear equations (12. 2t at a given Rayleigh 
number to create the nonlinear axisymmetric solution (U, H), we integrated the non-axisymmetric 
linearised equations ( 12.61 ) to evolve (u,/i) from an arbitrary initial condition. To integrate J2.6I ), 
we used a timestep of At = 10~ 4 and a spatial resolution of N r =47,N Z = 29 for each azimuthal 
mode. To construct the Krylov space ( 12.91 ) and approximate eigenpairs, we used K = 10 vectors, 
a time interval of T = lOOAf = 10~ 2 , and an acceptance criterion of 10 -5 . 

2.5. Complex eigenvectors and their representations 

The linear problem ( 12.61 ) for perturbations (u,h) about an axisymmetric convective state (U,H) 
can be divided into decoupled subproblems, each corresponding to a single azimuthal wavenum- 
ber m. The problem for wavenumber m can in turn be divided into two identical decoupled 
subproblems, corresponding to fields of the form 

u r (r,z) cos (mQ), u%(r,z) sin(m9), u z (r,z)cos(mQ), h(r,z) cos(mG), (2.12a) 

and 

u r (r,z) sin(m8), we(r,z)cos(w8), u z (r,z) sin(m9), «(r,z)sin(m8). (2.12b) 

For simplicity, we will represent each of these types of vector fields by its temperature component 
h(r,z) and leave the dependence on and on t to be written explicitly. We may write the linear 
evolution problem (12.71 ) restricted to fields with trigonometric dependence on mQ such as d2.12al >- 
( 12.12/?! ) as 

d,h = L m h. (2.13) 

A real eigenvalue breaking azimuthal symmetry in an 0(2) symmetric situation is associated 
with a two-dimensional eigenspace, consisting of linear combinations of vectors of type d2.12q| ) 
and §2.\2b\ . Since 

ah(r,z) cos(m8) + $h(r,z) sin(m8) =C^(r,z)cos(m(8-8 )), (2.14a) 

where 

C = yV + P 2 , mQ Q = atan(p/a), (2.14i) 
all real eigenvectors have m nodal lines and reflection symmetry about some Qq. If we take 
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C °c y/Ra — Ra C 2 and add (12. 14at to the basic axisymmetric state, we obtain the "circle" of steady 
states resulting from a circle pitchfork mentioned in § 12.21 

A complex eigenvalue in the 0(2) symmetric situation is associated with a four-dimensional 
eigenspace. Within each eigenvector class ( 12.12ab and (12.12Z?b . the eigenspace is two-dimensional, 
spanned by two linearly independent eigenvectors h and h , which are transformed by L„, as 



h R \ _( ii -co \ / h R 



h 1 I \ co n J \ h 1 



(2.15) 



In ( I2.15l l, h can be replaced by any linear combination of h R and h 1 , but once h R is selected, the 
choice of h 1 follows from ( 12.15b . Although the components of equation ( 12.15b are the real and 
imaginary parts of the complex equation 

L m (h R + ih 1 ) = (ji+ i(0)(ti R + ih 1 ), (2.16) 

the customary designation of h and h 1 as the real and the imaginary part of the eigenvector is 
arbitrary, as reflected by the fact that an eigenvector can be multiplied by any complex number. 

To form eigenvectors of the full cylindrical problem belonging to the four-dimensional eigen- 
space, each of h and h is multiplied by a trigonometric function. This yields as a basis for the 
four-dimensional eigenspace: 

h R (r,z) cos(w6), (2.17a) 
h 1 \r,z) cos(m8), (2.17b) 
h R (r,z)sm(tnd), (2.17c) 
V{r,z) sin(m8). (2.\ld) 

One choice for a complex eigenvector pair is (12.17ab -( f2.17frb . since 



h R cos(md) \ _ ( n -co \ f h R cos(mQ) 
h'cosimB) J ~ \ CO /J J \ /Vcos(m6) 



(2.18) 



More generally, the trigonometric dependence can be taken as in d2.14ab . with the same trigono- 
metric dependence for each of h and h 1 , to form a complex conjugate eigenvector pair each of 
whose members has m nodal lines and m axes of reflection symmetry, including 8 = 0o. The 
evolution in time under ( I2.13l l for a field with an initial condition of this form is 

h{r,e,z,t)=ae fa [/i K (r,z)cos(cof)-^ / (r,z)sin(cof)]cos(m(e-9o)). (2.19) 

The subspace of fields with azimuthal dependence cos(m(0 — Go)) is invariant under linearised 
time evolution. (There also exists an invariant subspace under the nonlinear time evolution, which 
includes harmonics cos(km(Q — 8o)), with the same m axes of reflection symmetry.) If we take 
/j = and a <* ^/Ra — Ra c 2 in ( 12.191 ). and add this to the basic axisymmetric solution, then we 
obtain to first order the standing wave solution mentioned in § 12.21 

Any combination of d2.17ab - (l2.17a1 > is also a member of a complex eigenvector pair. The 
calculation 

uh R (r,z)cos(mQ) + pA 7 (r,z) sin(m8) 
aWfaz) cos(»?8) — $h R (r,z)sm(mQ) 

/j —CO A / a/i fi (r,z) cos (w8) + p/V(r,z) sin (m8) 
CO fi J \ a,h I (r,z)cos(mQ) — $h R (r,z)sin(mQ) 

when compared with ( 12.15b . shows that the two components of the vector in ( 12.20b form a com- 
plex conjugate pair of eigenvectors for the full cylindrical problem, as in ( 12.15b . Because h R (r,z) 
and h J (r,z) have different functional forms in (r,z), these vectors, unlike those of (I2.14ab . cannot 
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be combined into a single trigonometric function. Neither of the two components of (12.20b has 
nodal lines or reflection symmetry about any axis if both a and p are non-zero. The evolution in 
time under ( 12.131 ) for a field whose initial condition is the first component of ( 12.201 i is 

h(r,Q,z,t) = e A "[/i^(r,z)(acos(m9)cos(cof) - psin(ra0) sin(cof )) 

+h I (r, z) (a cos (m6) sin (cof ) ) + p sin(ra0) cos (cof ) )] . (2.21) 

If p = ±a, then ( I2.2U becomes 

h(r,B,z,t) ^e f "a[h R {r 7 z)cos(mQ±m)+h I (r,z)sm{mQ±(S)t)}, (2.22) 

where f or 8 may be replaced by (t — to) or (0 — 0o). If we take /j = and a \/Ra — Ra c 2 in 
( 12.221 ) and add the basic axisymmetric solution, then we obtain, to first order, the expression for 
clockwise (mQ + cof) or counterclockwise (mQ — cof ) travelling waves mentioned in § 12.21 



2.6. Amplitude equations and normal form 

The linearised evolution treated in the previous section permits any combinations of G.llaj - 
il.lldj . T he mathemat i cal analysis of Hopf bifurcation in the presence of 0(2) symmetry carried 
out bv e g.lBaiaildl982h:lGolubitskv & Stewartldl985h : lKnoblochl(ll986l) : lvan Gils & Mallet-Paretl 



out by e.g 
([l986);K 



Kuznetsov ( 1998) describes the effect of including generic nonlinear terms compatible 
with the symmetries. Following the formulation of these authors, we decompose the field into a 
sum of clockwise and counterclockwise travelling waves with complex amplitudes = p-e'^~ 
and = p + e'^ + , respectively. The four variables p±,0± form another description of the four- 
dimensional space described in the previous section. The nonlinear evolution of £± near the 
bifurcation can be described by the following amplitude equations or normal form: 

c+ = (A<+/co+ fl |c-i 2 +Hii;+i 2 +ic-i 2 )K+, ( 2 - 23 «) 

t = (^ + /co + fl |C + | 2 + KIC + | 2 + IC-| 2 ))C- (2.236) 

We use the normal form to interpret the results of our full numerical simulations. 
Separating d2.23t into equations for real amplitudes p± and phases §± leads to 

p + = (fi + ehpl+b r (p\ + pi))p+, (2.24a) 

p_ = (/j + a r p 2 + +b r (p 2 + + p 2 _))p-, (2.24b) 

<j)+ =(ti + ajp 2 _+bi(p 2 + + p 2 ), (2.24c) 

§-=-(£>-aip 2 + -bi(p 2 + + p 2 _). (2.24d) 

Periodic solutions to d2.24t must be either standing or travelling waves. Solutions to d2.24t and 
their properties are given in Table Q] This table shows that both standing and travelling wave 
solutions exist for /j > if b r and a r + 2b r are both negative. A positive growth rate from a 
solution indicates instability. Thus, the stability of the solutions depends on the sign of a r : if 
a r > 0, then standing waves are stable and travelling waves unstable, and vice versa for a r < 0. 
Figure |2] shows phase portraits for the amplitudes (p+,p_), for the cases in which all three 
branches co-exist and either the standing or the travelling waves are stable. 



3. Results 



3.1. Conductive state 

Figure [3] shows the linear stability limits of the conductive state to perturbations with azimu- 
thal wavenumbers m = 0, 1, and 2 (Bor onska & Boro riski 2001). These results, obtained with 



the linearised version of our code, agree very closely with those presented by Wan schura et al. 
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solution growth rates frequencies 



Basic state p + = p_ = /j, /j 

Counterclockwise wave p + = \J~^-, P- = 

Clockwise wave p_ = \ / ^ , p+ = 



-f 



Standing wave p+ = p_ = v 

TABLE 1 . Solutions to l |2.24t and their properties 







a,. 




2a r 






o >o< — o >•< — 

TW+ TW + 

FlGURE 2. Phase diagram illustrating stability of standing waves (left) or travelling waves (right). The 
origin is the basic state and the axes represent amplitudes of counterclockwise and clockwise travelling 
waves p+ and p_ . Standing waves can be constructed as an equal superposition of the two. 
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FIGURE 3. Linear stability of the conductive state 



(1996). Note that in the range 0.9 < T < 1.57, the primary instability is axisymmetric. Imme- 
diately below and above this range of aspect ratio, the first instability is to an eigenvector with 
azimuthal wavenumber m=l. Instability of the conductive state is independent of Pr. However, 
the resulting nonlinear states and their stability depend on Pr; in the remainder of the study we 
fix Pr = 1. 




FIGURE 4. Temperature contours for axisymmetric solutions at T = 1.47 and Ra = 1950 with upward (left) 
and downward (right) flow at the centre. Solid (dashed) curves correspond to positive (negative) values, 
here and in subsequent visualisations. 



3.2. Steady axisymmetric state 

W e reproduced the primary flow for r = 1.47 and Ra = 1950, parameters for which, according 
tolWanschura et al.\ dl996l) and figure [3] the conductive state is unstable only to axisymmetric 
perturbations. In a fully three-dimensional simulation, starting the evolution from an arbitrary 
non-axisymmetric perturbation about the conductive state, we obtained a flow consisting of one 
toroidal roll. While axisymmetric, this flow breaks the reflection symmetry in z and thus two such 
states exist, with either upflow or downflow at the centre; these are illustrated in figure [4] We 
used the state with downflow at the centre as the initial condition for higher Rayleigh numbers. 
According to the calculations of Wanschura et al. ( 1996), the axisymmetric state first bifurcates 



towards a flow with azimuthal wavenumber m = 3 for 1.45 < T < 1.53 and with wavenumber 
m = 4 for 1.53 < r < 1.57. The critical Rayleigh numbers Ra c i at which this loss of stability 
occurs are given in table [2] 

3.3. Eigenvalues and eigenvectors 

Using the methods described in § 12.41 we integrated the evolution equations ( 12.61 > linearised about 
axisymmetric solutions for aspect ratios 1 .45 < T < 1 .57 and several different Rayleigh numbers. 
The leading eigenpairs calculated for Ra = 24000, F — 1.57 are given in Table [3] For these 
parameter values, the critical eigenvectors are (in order of decreasing growth rate): two conjugate 
pairs with azimuthal wavelengths m — 4 and m ~ 3, a real eigenvector with m = 1, and another 
conjugate pair with m = 5. 

Figure [5] represents the dependence of the leading eigenvalues on Rayleigh number for as- 
pect ratios T = 1 .47 and T = 1 .57, along with the azimuthal wavenumbers of the corresponding 
eigenvectors. Ra C 2 was calculated by determining the zero crossing of /j (Ra), the growth rate of 
the leading eigenvalue (that of largest real part), by linear interpolation. (Critical Rayleigh num- 
bers calculated by introducing perturbations into nonlinear simulations at various values of Ra, 
and fitting the initial evolution to an exponential to calculate growth or decay rates fj(Ra) gave 
similar results.) We then calculated co C 2 = Q>(Rac2), also by linear interpolation. The values we 
ob tained for two aspect rat ios T = 1 .47 and F = 1 .57, and the corresponding values published 



sp ect rat i 

bv lWanschura et al.\ ( 1996) are those given in Table [2] The critical wavenumbers are the same, 



and the errors in Ra C 2 and in co c 2 are less than 1%. In what follows, we will focus on the m = 3 
instability, since the m = 4 transition is similar; the aspect ratio is F = 1 .47 unless otherwise 
specified. 

We summarise here the differences between our numerical method and that of Wanschura 
et al. (1996). We linearised a timestepping code in order to, in effect, carry out the power 
method (supplemented by an Arnoldi decomposition) on the exponential exp(LAf) of the Ja- 
cobian. Wanschura et al. constructed the Jacobian matrix L and used inverse iteration to compute 
its eigenvalues. Our calculation was restricted to one of the two identical decoupled subprob- 
lems, corresponding to only one of the invariant subspaces of the form (I2.12al i or (I2.12fcb . As a 
result, the complex eigenfunctions we show in tableware all in the eigenspace corresponding to 
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r present study Wanschura et al. error 



Ra cl 24738 24928 0.76% 

1.47 a) e2 42.33 42.54 0.48% 

m c 2 3 3 



Ra c2 22849 23011 0.70% 

1.57 co e2 45.26 45.47 0.45% 

m c 2 4 4 

TABLE 2. The parameters of the oscillatory bifurcations found by linear analysis: critical Rayleigh 
numbers Ra c 2, critical frequencies 00^2 an d azimuthal wavenumbers of critical eigenvectors for two aspect 

ratios. 



, eigenvector visualisation . 

eigenvalue f ... wavenumber 

real part ± imaginary part 

0.86±46.3( 4 



10 



-10 



0.24^41.6/ OM&O gig 



3 2x10" 



-0.81 



1 6x 10 



-10 



-4.40-45.9/ Qi 



Cu f 



9 x 10" 07 



TABLE 3. For Ra = 24000, T = 1.57: eigenvalues, visualisation of corresponding eigenvectors, azimu- 
thal wavenumber and residual error. The visualised field is the temperature at the midplane; for complex 
conjugate eigenpairs the real and imaginary parts of the eigenvector are depicted. 



standing waves, with three axes of reflection symmetry. Basis vectors for the remainder of the 
four-dimensional eigenspace can be found by rotating the eigenvectors of table[3] i.e. multiplying 
by sin(m0) instead of cos(m8). Wanschura et al, in contrast, used the travelling wave form as an 
initial condition or invariant subspace, as discussed below. 

In figure [6] we show representative elements of the eigenspace associated with the m = 3 
complex eigenvector at Ra = 25000. Figures [6](a, b) show h R (r,z)cos(mQ) and h (r,z)cos(m8), 
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Ra 
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FIGURE 5. Leading eigenvalues as a function of Rayleigh number for aspect ratio T = 1.47: (a) real part, 
(b) imaginary part and for aspect ratio T = 1.57: (c) real part, (d) imaginary part. Vertical thin dashed line 
marks Ra c2 = 24738 for T = 1 .47 and Ra c2 = 22849 for T = 1 .57. 




FIGURE 6. Eigenvectors for T = 1.47, Ra = 25000 (temperature field contours at z = 0): (a) real part of the 
critical eigenvector; (b) imaginary part of the critical eigenvector; (c-g) superposition of the two fields via 
h R (r,z) cos(m(8 - G )) +h'(r,z) sin(m(6 + 9 )), with m% of (c) 0, (d) jc/4, (e) it/2, if ) 3jt/4, (g) 0.9271. 



while figures |6](c-g) are generated via 

C (h R (r, z) cos(m(6 - 6 )) + V(r, z) sin(m(6 + Go))) , (3.1) 

a form equivalent to (12. 2U after translation of 8 and of t. Clockwise travelling waves ensue 
for mQo — k/2 (c), counterclockwise travelling waves for m8o = (e), and standing waves at 
different temporal phases for mQo = ±Jt/4 (d,f). Thus, the angle mQo is similar to that used in 
figure |2] An eigenvector which corresponds to neither travelling nor standing waves is shown in 
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figure |6](,g). These are all depicted on the slice z = 0; when we plot the field of figure |6fc) at 
z = 0.3, we recover the form shown by Wanschura et al. We emphasise, however, that the other 
fields depicted in figure [6] are all equally valid eigenvectors. In particular, a nonlinear analysis, 
such as the simulations presented below, is required to determine whether the resulting nonlinear 
flow near onset is a travelling or a standing wave. 

3.4. Weakly unstable standing waves 

Above the critical Rayleigh number Ra C 2, a slightly perturbed axisymmetric state evolved in our 
simulations towards a three-dimensional time-dependent state, presented in figures|7l[8]and|9] Fi- 
gure|7]shows temperature contours on the midplane at six regularly spaced instants in time within 
one oscillation period. In contrast to the eigenvectors depicted previously, figure [7] displays full 
nonlinear temperature fields, which are dominated by a large axisymmetric component. There 
are six pulsing extrema, engendering oscillation between two triangular structures of opposite 
phases (figures|7]fl and[7]<i). At each instant, the flow is invariant under rotation in 8 by 2n/3. In 
addition, this flow is also symmetric with respect to three different axes of reflection. Figure [8] 
shows contours of azimuthal velocity at the same times as figure [7] Figure [9] shows the tempe- 
rature dependence on the angle 9 for fixed radius and height at different times. Six fixed nodes 
identify this state as a standing wave with azimuthal wavelength 2jt/3. 

The standing wave state persists for such a long time that it might seem stable. However, a 
small reflection-symmetry breaking imperfection develops that eventually leads to the transition 
to travelling waves. Figure [10] shows the temperature dependence on the angle for the same 
parameters as figure|9] but at a later time. The breaking of reflection symmetry can be observed 
when the amplitude of the standing wave is small. The standing waves can be stabilised by 
imposing reflection symmetry. When we did this, above a threshold Ra C T, sa 27 000, we discovered 
a new (unstable) standing-wave solution, displayed in figure QT|for R a = 30000. 

In order to study the transition from standing to travelling waves, we monitored the growth 
of antisymmetric components. When the standing wave is still dominant, the amplitude of the 
antisymmetric components behaves in time like (A cos (Of +B) exp (jigw-ytw t), where fJ S w^tw is the 
growth rate from standing waves to travelling waves. The growth rate fJ sw ^ t w, shown on figure[T2l 
as a function of Ra, is about two thirds of fJo^i, the growth rate from the axisymmetric state to an 
m = 3 flow (denoted in the previous sections by fj). The observed lifetime of the standing waves 
decreases as the Rayleigh number is increased, since the growth rate fJ S w^tw increases. 

3.5. Stable travelling waves 

After the pattern has evolved sufficiently from the standing wave state, the fixed antinodes 
abruptly begin to rotate about the cylinder axis. The six pulsing spots change into three rotat- 
ing spots, as the standing waves become travelling waves with the same azimuthal wavelength. 
Figure Q~3] [14] and [15] depict temperature profiles and contours of the temperature and the azi- 
muthal velocity of the travelling waves at different times. The travelling waves, like the standing 
waves, have three-fold rotational symmetry, but do not have reflection symmetry. 

Travelling waves are the final state of the time evolution. The reason for which we obtained 
standing waves before travelling waves in our simulations is that our initial conditions were 
reflection symmetric and our numerical procedures introduce antisym metric perturbat i ons at a 
low rate. (This is also seen in the simulations of thermocapillary flow by lLeypoldt et al . (2000).) 
When the Rayleigh number is decreased, travelling waves persist until Ra reaches Ra c i. 

We conducted simulations for several values of F in the range 1 .45 < T < 1 .53 and observed 
weakly unstable standing waves and stable travelling waves for all of them. We believe that the 
same scenario also occurs for 1.53 < T < 1.57, but with azimuthal wavenumber m = 4 instead 
of m = 3. 
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(a) (b) (c) (d) (e) (f) 

FIGURE 7. Standing waves at ifa = 26000: temperature contours on the midplane at t = 0, T/6, 2T/6, . . . 




(a) (fe) (c) W (e) (f) 



FIGURE 8. Standing waves at Ra = 26000: contours of azimuthal velocity on the midplane at t = 0, T/6, 

XT '/6,... 




FIGURE 9. Standing waves at = 26000: temperature versus 6 at (r,z) = (0.7,0.3) at five successive 

times. 




3.6. Amplitudes and frequencies 
We calculated the energy E of both types of waves by first defining a norm whose square is 



Ra V Pr Ra 
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(a) (b) (c) {d) (e) (f) 

FIGURE 1 1 . Oscillatory solution obtained at Ra = 30000 by imposing reflection symmetry: temperature 
contours on the midplane at / = 0, T/6, 2T/6, . . . 
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FIGURE 12. Growth rates as a function of Rayleigh number. Solid line: growth rate /jo^t, of m = 3 eigen- 
vector (either standing or travelling waves) from the axisymmetric solution (from linear evolution). Squares: 
growth rate /jsw^tw of travelling waves from standing waves (from nonlinear simulation) with linear fit as 
dashed line. 



where (,) denotes spatial integration; (13. 2\ is one of many possible choices for this system. We 
then simulated the nonlinear evolution equations and calculated (u,h) as the difference between 
the three-dimensional and the axisymmetric solution. We define E to be the integral of ( 13.2b over 
one oscillation period. 

The energies E sw , Ef W and frequencies C0 JM ;, Q&tw as function of Ra are shown in figure [161 
The energies and frequencies for the two types of waves are quite close. The frequency COo^3 
obtained from linear stability analysis is also reproduced from figure |5](£>) for comparison. For 
both types of waves, the frequencies near the threshold are close to the Hopf frequency and the 
energy satisfies E °< (Ra —Ra C 2). These are hallmarks of a supercritical Hopf bifurcation. 

3.7. Normal form coefficients 

Using the growth rates, amplitudes and frequencies of the standing and travelling waves that we 
have presented in sections 13.31 and 13.61 it is possible to calculate the coefficients of the normal 
form J2.24t for our particular case. The bifurcation parameter /j = /jo^3 and frequency CO = 0)0^3 
vary linearly with Ra —Ra C 2, while the other coefficients a r , b r , a,, bj are constants. 
From the data in figures |5](a,Z?), we extract the fits 

Ra — Ra c 2 



From the data in figure [ 



= 14.98 

000^3=42.33+21.21 
'lwe extract the fits 



Ra C 2 
Ra — Ra c 2 
Ra c 2 



E tw =Af w = p z + = -f =0.2037 



Ra — Ra c 2 
Ra c 2 



(3.3a) 
(3.3b) 

(3.4a) 
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(a) (b) (c) {d) (e) (f) 

FIGURE 14. Counterclockwise travelling wave at Ra = 26000: temperature contours on the midplane at 

t = 0, r/6,2T/6, ... 




FIGURE 15. Counterclockwise travelling wave at Ra = 26000: contours of azimuthal velocity on the 

midplane at t = 0, T/6, 2T/6, . . . 



o o o — u Ra — RcifO 

E sw = A sw = p 2 + + p 2 _ = 2—^— = 0.13 — 2, (3.46) 

a r + Lb r Ra c 2 

«W = ©o^i -t-H = 42-33 + 16.26 Ra ~ Ra<:2 ^ ( 3-4c ) 
6 r /?a C 2 

a,- + 26; Ra — Ra r ? 
co m . = co ^3 - i^tTT" P = 42.33 + 17.29 — 2. {3Ad) 

Equations ( 13.4b are used to determine the nonlinear coefficients as 

6 r =-73.5, (3.5a) 

a r = -83.6, (3.56) 

bi = -24.3, (3.5c) 

Oi = 11.7. (3.5*/) 

An additional equation is provided by the data in figure Q~2] showing the growth rate fi sw ^tw from 
standing to travelling waves: 

2a r Ra — Ra C 2 
Vsw^tw = — ——^=10.23 — . (3.6) 

a r + 2b r Ra C 2 
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FIGURE 16. Dependence of energy and frequency on Rayleigh number for standing and travelling waves. 
Vertical dashed line indicates the critical Rayleigh number Ra c 2 for onset of the waves. 



and provides a second determination of a r 



-78.8. 



(3.7) 



which differs by 6% from i3.5b\ . 



4. Conclusion 

We have used both nonlinear simulations and linear stability analysis to elucidate the be- 
haviour of Rayl eigh-Benard convection in the parameter region of 1.45 < T < 1.57, Pr = 1 
first studied by IWanschura et al.\ (I1996T) . In this regime, the primary axisymmetric convective 
state loses stability to an m = 3 perturbation via a Hopf bifurcation whose critical eigenspace 
is four-dimensional. We calculated representative eigenvectors and explained how these relate 
to those computed by Wanschura et al. The bifurcation scenario guarantees that branches of 
standing waves and of travelling waves are created at the bifurcation, but that at most one of 
these branches is stable. Our nonlinear simulations showed a supercritical bifurcation leading to 
long-lived standing waves which were eventually succeeded by travelling waves, both as time 
progressed and as the Rayleigh number was increased. We explained this by showing that the 
rate of transition from standing waves to travelling waves, while positive, is nevertheless small. 
In the absence of long-time integration and of these analyses, it would be easy to conclude that 
the standing waves were stable. This underlines the importance of calculating growth rates, in 
addition to carrying out nonlinear simulations, and of using established bifurcation scenarios to 
interpret physical phenomena. 

The numerical and theoretical techniques we have used can be generally applied to study 
transitions in hydrodynamic problems. Our main tool was direct numerical simulation of the 
governing Boussinesq equations using a pseudo-spectral semi-implicit timestepping code. We 
complemented this approach with several other techniques. To carry out stability analysis, we 
first linearised the code. This requires very little modification of the existing code, but yields 
results which are far more precise and robust than restricting integration to the time interval du- 
ring which perturbations to the basic state are small. Integrating the linearised equations is, in 
effect, an implementation of the power method for finding the fastest growing eigenvalues and 
corresponding eigenvectors. Eigenvectors with different azimuthal wavenumbers can be found 
simultaneously, since the linearised evolution of each Fourier mode is independent of the others. 
For a single wavenumber, this use of the power method is rendered more accurate and more 
general by postprocessing the results of linearised time integration with the Arnoldi decomposi- 
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tion to extract several, possibly complex, eigenvectors. We also interpreted our results in light of 
known results concerning axisymmetry-breaking Hopf bifurcations in systems with 0(2) sym- 
metry. This framework allows us to generate the four-dimensional eigenspace by combining 
eigenvectors with different symmetries. Traditionally, eigenvectors corresponding to clockwise 
and counterclockwise travelling waves are combined to form standing waves; we used a com- 
plementary, but equivalent, approach of combining standing waves of different spatial phases to 
form travelling waves. Finally, we interpreted our results in terms of the four ordinary differen- 
tial equations comprising the normal form for Hopf bifurcations in systems with 0(2) symmetry. 
Using our nonlinear simulations of the governing Boussinesq equations, we were able to calcu- 
late the various coefficients in the normal form equations. 

We have not sought to determine the limits of the range of this phenomenon, in aspect ratio 
and Prandtl number. As these ranges were given by Wanschura et al. only for Pr = 1, a future 
direction would be to determine the whole zone in the parameter space where the Hopf bifur- 
cat ion occurs . It wo uld be interesting also to examine more closely the pulsing pattern found 



bv lHof et ai.l (119991) at Ra = 33000, T = 2, Pr = 6.7, in order to determine whether this state, 
evolving from axisymmetric flow, is the result of a bifurcation similar to that described in the 
present paper. 

The computations were performed on the NEC SX5 of the IDRIS (Institut du Developpement 
et des Ressources en Informatique Scientifique) supercomputer center of the CNRS (Centre Na- 
tional pour la Recherche Scientifique) under project 1119. 
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